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ABSTRACT 
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CN ■ The direct detection of gravitational waves is a major goal of current astro- 

physics. We provide details of a new method for detecting a stochastic back- 
ground of gravitational waves using pulsar timing data. Our results show that 
regular timing observations of 40 pulsars each with a timing accuracy of 100 ns 
will be able to make a direct detection of the predicted stochastic background 
from coalescing black holes within five years. With an improved pre-whitening 
algorithm, or if the background is at the upper end of the predicted range, a 
*0 ■ significant detection should be possible with only 20 pulsars. 

o 
in 
o 

6 

Analysis of pulsar pulse time-of-arrival (TOA) data shows that pulsars, especially mil- 
lisecond pulsars (MSPs), are very stable clocks. Measurement of timing residuals, that 
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Introduction 



is, the differences between observed and predicted TOAs, enables the direct detection of 
gravitational waves (GWs) (Estabrook & Wahlquist 1975; Sazhin 1978; Detweiler 1979). 
The fluctuating TOAs induced by a GW will be correlated between widely-spaced pulsars. 
Hellings & Downs (1983) attempted to detect this correlation by cross-correlating the time 
derivative of the timing residuals for multiple pulsars. In our work, we have developed a 
similar cross-correlation technique and have, for the first time, a fully analyzed method for 
combining multiple pulsar observations in order to make an unambiguous detection of a GW 
background. We emphasize that, in contrast to Hellings & Downs (1983), our method is 
based entirely on the measured residuals. 
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Only the effects of a stochastic background of GWs are considered. Astrophysical sources 
of such a background include cosmological processes (e.g. Maggiore 2000) and coalescing 
massive black hole binary systems (Jaffe & Backer 2003; Wyithe & Loeb 2003; Enoki et 
al. 2004). We show that a direct detection of a stochastic GW background is possible 
using pulsar timing observations and that the significance of the detection depends upon the 
number of pulsars observed, the root-mean-square (RMS) timing noise achieved, the number 
of observations, and the power spectrum of the measured timing residuals. The results are 
applied to the case of the Parkes pulsar timing array (PPTA 1 ). 

In the next section, the analysis technique is described. In §3 the significance of detecting 
a given stochastic background using this method is estimated. The effects of pre-filtering 
the residual time series are also discussed. The results are summarized in §4. 



2. Detection Technique 

As a first step, the power spectra of the pulsar timing residuals are analyzed. If they all 
show a very red power-law spectrum, the residuals may be dominated by a GW background. 
However, such red spectra can also be due to period noise intrinsic to the pulsar, uncorrected 
interstellar delays, inaccuracies in the Solar-System ephemeris, or variations in terrestrial 
time standards (e.g. Foster & Backer 1990). A GW background produces a unique signature 
in the timing residuals which can only be confirmed by observing correlated signals between 
multiple pulsars widely distributed on the sky. 

The presence of a stochastic GW background will cause the pulse TOAs to fluctuate 
randomly, but these fluctuations will be correlated between different pulsars. In order to de- 
tect the presence of a GW background, one needs to first calculate the correlation coefficient 
between the observed timing residuals of each pair of observed pulsars: 

1 

where R(ti, k) is the time series of N pulsar residuals sampled regularly in time, k\ and k 2 
are the directions to the two pulsars, and cos(#) = k\ ■ k,2- It will be assumed that R has zero 
mean and that each pulsar pair has a unique angular separation. r{9) is written only as a 
function of the angular separation since the GW background is expected to be isotropic. In 
the presence of an isotropic GW background, the ensemble- averaged value of r{6) is given 



1 http: / / www.atnf.csiro.au/research/pulsar /psrtime 
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by 2 : 

(r(0)> = aim (2) 
C(0) = \x\og{x)- X -+ l -+ l -5{x) (3) 

where x — (1 — cos(0))/2, cr 9 is the RMS of the timing residuals induced by the stochastic 
GW background, and S(x) equals 1 for x = and otherwise. The detection technique 
proposed here simply looks for the presence of the function ((9) in the measured correlation 
coefficients r{9). 

Since one cannot perform the ensemble average in practice, the measured statistic, r(9), 
will be of the form r(6) = (r(0)) + Ar(9), where Ar(9) is a "noise term". Since r(9) is 
calculated by summing over a large (> 20) number of data points, Ar{9) will be a Gaussian 
random variable for practical purposes. The optimal way to detect the presence of a known 
functional form within random data is to calculate the correlation between the data and the 
known function. Hence, to detect the presence of the GW background one needs to calculate 

iEivV(^-r)(C(^)-C) 

p = ( 4 ) 

where 9i is the angle between the ith pair of pulsars and N p is the number of distinct pairs 
of pulsars, f and ( indicate the mean values over all pairs of pulsars and of and <r 2 are the 
variances of r and ( respectively. For M pulsars, N p = M(M — l)/2. 

From the definition of r(9) and Eqn. 4, one can show that the expected value of p is 
approximately: 



(5) 



Np-l 

< = ttE«^)-^))) 2 )- ( 6 ) 



p i=0 



For the case where there is no correlation in the data, the statistics of p will be Gaussian 
with zero mean and variance given by a 2 p = 1/N P = 2/(M 2 — M). Hence, the significance of a 
measured value of p may be defined as S = p/cr p . The probability of measuring a correlation 
greater than or equal to p when no actual correlation is present is given by erf(S'/v / 2)/2. 



2 For an outline of the calculation of £ see Hellings & Downs (1983). 
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3. Estimating the Detection Significance 

In order to estimate the expected detection significance, S, one needs to estimate a g 
and oat- It is assumed that the timing residuals, R(t,k), are stationary Gaussian random 
variables that are sampled at regular intervals denoted by At. It is also assumed that terms 
proportional to t and t 2 (i.e., the period and period-derivative terms) have been subtracted 
from R(t, k). 

The space-time fluctuations induced by a stochastic GW background are described by 
a quantity known as the characteristic strain spectrum denoted by h c (e.g. Maggiore 2000). 
Models of the GW background propose a power-law dependence between h c and the GW 
frequency, /: h c (f) = Af a (Jaffe & Backer 2003; Wyithe & Loeb 2003; Maggiore 2000; Enoki 
et al. 2004). Using this form of the characteristic strain spectrum, the power spectrum of 
the induced residuals is given by -Pr(/) = (\R(f)\ 2 ) = ^j/ 2 " -3 , where R(f) is the Fourier 
transform of R(t). Given Pr(/), the total RMS fluctuation induced by the stochastic GW 
background is given by 

fh 

" PrUW (7) 
A 2 

~ 2tt 2 (2 - 2a) [Jl Jh ] K ) 

where f\ is the lowest detectable frequency given by 1/T and fh is the highest detectable 
frequency typically given by 1/2 At. T is the total time span of the data set. Since a < 
for backgrounds of interest (Maggiore 2000), the term containing fh is negligible. 

Estimating o /\ r is slightly more complicated. To take into account the effects of subtract- 
ing linear and quadratic terms from the residuals, a semi-analytic approach was adopted. As 
outlined below, an estimate for a^ r is made analytically but with one free parameter (3. For 
a given value of f3, S is calculated as a function of A for a given set of pulsars and timing 
parameters. S(A) is compared to Monte-Carlo simulations in order to determine the correct 
value of p. This showed that the value of (3 is insensitive to the values a, N, M, a g and the 
RMS residual noise level. 

Using Equation 1 together with the assumption that R(t, k) is a Gaussian random 
variable, one can show that 

~ N-l N-l 

^Ar ~ JpJ2J2 Ci ^l) C ij(k2), (9) 
i=0 j=0 

where Cij(k) = (R(t + iAt, k)R(t + jAt, k)). The bar above Equation 9 represents an average 
over all pairs of pulsars. As the autocorrelation function and the power spectrum are Fourier 
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transforms of one another, one can estimate o\ r from the expected power spectrum of the 
residuals. The statistics of the residuals are assumed to be stationary so that <%•(&;) depends 
only oni- j. The expected discrete power spectrum of R(t, k), which includes both a GW 
component and a white noise component, is given by 

10 for % = 

P g (i) is the discrete power spectrum of the GW-induced timing residuals, i is the discrete 
frequency bin number corresponding to frequency i/T, <7 n (k) is the RMS value of the residual 
fluctuations caused by all non-GW sources for the pulsar in the k direction. It is assumed that 
all noise sources have a flat spectrum. This assumption is consistent with most observations 
of MSPs. P g {i) is given by 

^2y2-2a 

Pti) = - — — -m(i) (11) 

flW (27r) 2 (2-2a) w K J 

where 

m = for % — 

m = p2a-2 _ ( 15 )2a-2 for ! = 1 

m = (i- 0.5) 2a - 2 - (i + 0.5) 2q " 2 for i > 1. 

Effectively, (3 is the lowest frequency used to calculate the correlation function c^. Monte- 
Carlo simulations show that (3 = 0.97. 

For the case where all pulsars have the same noise level, the detection significance 
becomes 



S = 



M(M - l)/2 



1 + 



where x — ^r)v Silo 1 SjLo 1 C 9ij ' anc ^ c ^ ^ s ^ e correlation function for the GW-induced 
component of the timing residuals, x is a measure of the "whiteness" of the residuals. 

The solid curve in Figure 1 panel A) plots the detection significance versus power-law 
amplitude for a = —2/3, the expected value for a background generated by an ensemble of 
super- massive black hole binaries (Jaffe & Backer 2003). This spectral index together with 
the removal of the linear and quadratic terms from R effectively makes x — 0.6iV. The 
parameters are set as follows: iV = 250, M = 20, a n = 100 ns and T = 5 years. These values 
are the target values for the PPTA (Hobbs 2004). Note that the significance saturates for 
high values of A. This effect can easily be seen in Equation 12 since all terms of the form 
cr n /cr g go to zero as a g gets very large. This saturation is due to the "self-noise" associated 
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with the stochastic nature of the background and its asymptotic value is independent of a n . 
The roll-off at low values of A occurs at a g — a n . 

Since the power spectrum of the GW-induced timing residuals will be dominated by 
low frequencies, one can apply a low-pass filter to each of the residual time series before 
correlating. This is similar to fitting a low-order polynomial to the data and then correlating 
the resulting fits. To estimate the significance for this technique, one evaluates a and 
o\ r using Equations 8 and 10 but with a high frequency cut-off fh c . For purposes of this 
discussion, fhc was set to 4/T. The dashed line in Figure 1 panel A) shows the effect of using 
a low-pass filter on the residuals. All the other parameters are the same as those used to 
generate the solid line. Low-pass filtering effectively reduces a n while keeping a g relatively 
unchanged. It also has the effect of increasing x/N when P g is a red power-law spectrum. 
Hence, low-pass filtering will not increase the maximum attainable significance, but it will 
lower the value of a g where the roll-off starts to occur. 

We next try to increase the maximum achievable significance. This method involves 
both low-pass filtering and a technique called "whitening" . When correlating two time series 
that each have a steep power-law spectrum, an optimal signal-noise ratio is obtained if filters 
are applied to give each time series a flat spectrum before correlation. This will act to reduce 
X in Equation 12. In practice, starting from the lowest non-zero frequency bin, we give each 
Fourier component with significant power equal amplitude and set higher components to 
zero. In this way, we are correlating only that part of the signal which has a high signal-to- 
noise ratio and adjusting the power spectrum to optimize the measurement of the correlation 
function. 

Pd and a g need to be calculated in order to estimate S using the whitening method. 
After whitening P^{i,k) = 2cr ( i(k) 2 /N, where <Jd{k) is the RMS of the residual data from 
the kth pulsar. The whitening also affects a g . In the general case where the pulsars have 
different noise levels, a g will depend on the pulsar. The expression for p then becomes: 

p^ V ) = (13) 



with cr g (0) 2 given by 

a a( d ) 2 = i- T °d{k\)a d (k 2 ) 



2 



mi PS, h) £ P a (i)/P*(i, fci) ) (14) 

i=0 J \ i=0 



N 

where N max is the largest frequency bin used based on the criterion discussed above. The 
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solid line in Figure 1 panel B) plots the significance using the whitening versus A. The same 
parameters were used for this case as in the previous cases. 

The above discussion assumes that the noise levels were the same for all pulsars. Next, 
the case where the pulsars have different noise levels will be considered. All curves in Figure 1 
panel B) were generated using the whitening technique. Unless specified, 250 observations 
were taken on each pulsar over 5 years. The dashed line corresponds to 20 pulsars, 10 
with a n = 100 ns and 10 with a n = 500 ns. The dashed-dot line has 10 pulsars each with 
a n = 100 ns and 500 observations. The dashed-triple-dot line has 20 pulsars with a n — 100 ns 
and 500 observations over ten years. 

When given a choice between observing a large sample of pulsars with different noise 
levels and observing only those pulsars with the lowest noise levels but for a longer time, 
the above curves demonstrate that one should actually observe the larger sample of pulsars. 
This is not a general statement, but rather it depends on the level of the GW background 
and the noise level. However, the levels chosen above are relevant to the PPTA (Jaffe & 
Backer 2003; Wyithe & Loeb 2003; Hobbs 2004). Note that for large M, the significance 
scales as M. Hence, doubling the number of pulsars will double the expected significance. 



4. Summary 

The main goal of this work is to determine the effectiveness of an array of pulsars, such 
as the PPTA, for detecting a stochastic background of GWs. Using a simple correlation 
technique, the detection significance was calculated given the number of pulsars, the location 
of each pulsar, the TOA precision, the number of observations, the total time span of the 
data, and the amplitude and power-law index of the GW background. For the case where 
all pulsars have the same white-noise spectrum, Equation 12, may be used to calculate the 
detection significance. For the case of the PPTA, it was found that the maximum achievable 
significance will be about 3 for a background with spectral index a = —2/3 and A ~ 10~ 15 
which is the expected level of the GW background from an ensemble of super-massive binary 
black holes in galaxies (Jaffe & Backer 2003; Wyithe & Loeb 2003; Enoki et al. 2004). Note 
that lowering the RMS noise level will only decrease the minimum detectable value of A and 
not increase the maximum attainable significance. 

Low-pass filtering the timing residuals, or equivalently, fitting low-order polynomials 
(i.e. cubic terms) to the residuals and correlating the coefficients, does not increase the 
maximum attainable significance. The significance level is increased by pre-whitening of the 
timing residuals. Using whitening, it is estimated that the PPTA could obtain a detection 
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significance greater then 4 for A > 3 x 10~ 15 yr~ 2 / 3 provided that efficient whitening filters 
can be designed and implemented. This is an area of further study and will be addressed in a 
later paper. With the same qualifiers, increasing the total time span of the PPTA to 10 years 
would yield a significance greater then 4 for A > 10~ 15 yr~ 2 / 3 . Since the significance scales as 
the number of pulsars, doubling that number will double the expected significance. Hence, 
using the simple correlation technique described here without any pre-filtering, a stochastic 
background with A > lCT 15 yr~ 2 / 3 will be detectable at a significance of about 5.5 using 40 
pulsars observed 250 times over 5 years and each having 100 ns timing precision. 

Part of this research was carried out at the Jet Propulsion Laboratory, California Insti- 
tute of Technology, under a contract with the National Aeronautics and Space Administra- 
tion and funded through the internal Research and Technology Development program. The 
authors wish to thank Russell Edwards for useful discussions. 
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Fig. I. — The detection significance, S, versus the logarithm of the amplitude A of the 
characteristic strain amplitude h c (f). The strain spectral index a = —2/3, corresponding 
to an astrophysical background of GWs generated by super-massive binary black holes. 
The vertical lines bound the values of A expected by models of the GW background (Jaffe 
& Backer 2003; Wyithe & Loeb 2003; Enoki et al. 2004). In panel A), the curves were 
calculated with 20 pulsars each with RMS residual noise fluctuations of 100 ns. The solid 
line corresponds to the simple correlation technique. The dashed line includes the effect of 
low-pass filtering. Panel B) shows the effects of the whitening technique. The solid line was 
calculated with the same parameters as in A). The remaining curves were generated using 
different noise levels and number of pulsars. See text for further details. 



